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Abstract. A theory for the magnetization of ferromagnetic films is formulated 
within the framework of many-body Green's function theory which considers all com- 
ponents of the magnetization. The model Hamiltonian includes a Heisenberg term, 
an external magnetic field, a second- and fourth-order uniaxial single-ion anisotropy, 
and the magnetic dipole-dipole coupling. The single-ion anisotropy terms can be 
treated exactly by introducing higher-order Green's functions and subsequently tak- 
ing advantage of relations between products of spin operators which leads to an 
automatic closure of the hierarchy of the equations of motion for the Green's func- 
tions with respect to the anisotropy terms. This is an improvement on the method 
of our previous work, which treated the corresponding terms only approximately 
by decoupling them at the level of the lowest-order Green's functions. RPA-like 
approximations are used to decouple the exchange interaction terms in both the 
low-order and higher-order Green's functions. As a first numerical example we ap- 
ply the theory to a monolayer for spin S" = 1 in order to demonstrate the superiority 
of the present treatment of the anisotropy terms over the previous approximate 
decouplings. 

Keywords: Quantized spin model; Many-body Green's functions; thin ferromag- 
netic films. 
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1 Introduction 



There is increasing activity in experimental and theoretical investigations of the 
properties of thin magnetic films and multi-layers. Of particular interest is the 
magnetic reorientation transition which is measured as function of temperature and 
film thickness; for recent papers, see [0, ^ and references therein. 

The simplest theoretical approach to the treatment of thin ferromagnetic films 
in the Heisenberg model is the mean field theory (MFT), which can be applied 
either by diagonalization of a single-particle Hamiltonian[0 or by thermodynamic 
perturbation theory This approximation, however, completely neglects collec- 
tive excitations (spin waves), which are known to be much more important for the 
magnetic properties of 2D systems than for 3D bulk materials. In order to take the 
influence of collective excitations into account, one can turn to many-body Green's 
function theory (GFT), which allows reliable calculations over the entire range of 
temperature of interest: see, for example, Refs. [^, §,0], where the formalism in- 
cludes the magnetic reorientation. The application of Green's functions after a 
Holstein-Primakoff mapping to bosons, as applied in Ref. ||^, is valid only at low 
temperatures. Another method, which also can treat the magnetic reorientation for 
all temperatures, is the application of a Schwinger-Boson theory 0. Classical Monte 



Carlo calculations are also able to simulate the reorientation transition (see |T0[ and 
references therein). The temperature-dependent reorientation transition has also 
been investigated with a Hubbard model |jll| . 

In the present paper, we apply a Green's function theory to a Heisenberg Hamil- 
tonian plus anisotropy terms, a system previously treated at the level of the lowest- 
order Green's functions|^, ^,0. The approximate treatment of the single-ion anisotropy 
in the previous work is avoided here by extending the formalism to higher-order 
Green's functions and applying relations for products of spin operators, a procedure 
which leads to automatic closure of the hierarchy of equations of motion with respect 
to those terms stemming from the single-ion anisotropy. The exchange terms occur- 
ring in the higher-order Green's functions must, however, still be decoupled in an 
RPA-like fashion. This can be considered as an extension of the work of Devlin |1T2|, 



who has applied higher-order Green's functions to the description of bulk magnetic 
materials in one direction only. Our formulation applies to all spatial directions of 
a multi-layer system. We formulate the theory explicitly for a monolayer for spin 
S = 1 and give equations for an extension to the multi-layer case. It is straightfor- 
ward to see how the theory could be applied to higher spins. 

The paper is organized as follows: in Section 2, the previous theory^], ^ for thin 
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ferromagnetic films is generalized by introducing higher-order Green's functions, 
the model being explained in detail for a monolayer with S = 1. Section 3 gives 
the formal extension to the multi-layer case. In Section 4, numerical results for the 
S = 1 monolayer demonstrate the superiority of the exact treatment of the single-ion 
anisotropy term over the previously applied 0, ^] Anderson-Callen[0] decoupling. 
Section 5 contains a summary of the results. 

2 The Green's function formalism 

We investigate here a spin Hamiltonian, nearly the same as in Ref. consisting 
of an isotropic Heisenberg exchange interaction between nearest neighbour lattice 
sites, Jki, second- and fourth-order single-ion lattice anisotropics with strengths K2^k 
and K^^k respectively, a magnetic dipole coupling with strength qm, and an external 
magnetic field B = {B^, By, B'): 

^ = E JM{s,:st + sisn 

^ <kl> 
k k 

- E {^2^' St + \b+S, + B^Si) 

k 

+ \ E ^4{^USuSt + SlSt) - 3(S. ■ r«)(S, ■ r,,)) . (1) 

^ kl ''^kl 

Here the notation Sf = 5^ ± iSf and 5^ = B^ ± iB^ is introduced, where k and I 
are lattice site indices, and (kl) indicates summation over nearest neighbours only. 
Here, we add to the Hamiltonian in Ref. a fourth-order anisotropy term which we 
can treat exactly but for which we previously had no decoupling procedure available. 

Each layer is assumed to be ferromagnetically ordered: spins on each site in the 
same layer are parallel, whereas spins belonging to different layers need not be. Fur- 
thermore, the anisotropy strengths, coupling constants and magnetic moments are 
considered to be layer-dependent, so that inhomogeneous systems can be considered. 

To allow as general a formulation as possible (with an eye to a future study of 
the reorientation of the magnetization), we formulate the equations of motion for 
the Green's functions for all spatial directions: 
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G;^^^{uj) = {{Sr-Sj))^ (2) 
Gtfiu) = {{St;Sj))^. 

Instead of decoupling the corresponding equations of motion at this stage, as we 
did in our previous work[^, ^, we add equations for the next higher-order Green's 
functions: 



G^/'^M = {{i2S^-l)St;Sj))^ 
Gi/^^iu) = {{Sr{2S^-iy,SJ))^ 

Gj+'^M = {{StSt;Sj))^ (3) 

Gtrico) = mS!S!-2S{S + l));SJ))^. 

The particular form for the operators used in the definition of the Green's functions 
in Eqs. is dictated by expressions coming from the anisotropy terms. Terminating 
the hierarchy of the equations of motion at this level of the Green's functions results 
in an exact treatment of the anisotropy terms for spin S = 1, since the hierarchy 
for these terms breaks off at this stage, as will be shown. The exchange interaction 
terms, however, still have to be decoupled, which we do with RPA-like decouplings. 

For the treatment of arbitrary spin S, it is necessary to use 43(3 + 1) Green's 
functions in order to obtain an automatic break-off of the equations-of-motion hier- 
archy coming from the anisotropy terms. These are functions of the structure G'j'j'^ 
with a = {z)^{+)"^ and a = (— where, for a particular spin 3, all combina- 
tions of m and n satisfying (n + m) = 23 have to be taken into account. There occur 
no Green's functions having mixed + and — indices, because these can be reduced 
by the relation 3^3^ = 3{3 + 1) ^ - {3')'^. 

The equations of motion which determine the Green's functions are 

u;Gr;^ = A^f + (([0r,7^]_;^;)), (4) 
with the inhomogeneities 

A^^ = m,3j]^), (5) 

where Of are the operators occuring in the definition of the Green's functions, and 
(...) =Tr(...e-^^). 
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In the following, we treat a monolayer with S = 1 explicitly. In this case, a 
system of 8 equations of motion is necessary: 



k 

+(-^2,i + K4^^i)Glp^ — 2Ki^4^{G\j^ ^'^ — G^j^ ^'^), 
toGij^ = A^j^ + ^Jik{{S-S,^ — S^Si ; SJ)) — B^Gij'^ + B G^f^ 

k 

+ K4,)G,7'^ + 2K,4Gp''^"'^ - Gp 



k 

ujGlf = Alf + 2 X! Jik{{S^ - Sf. S^] SJ)) - -B Gp^ + -B^Gij'"^, 

k 



^G^^ = At;^^ -WjJ{{{6S!S^ - 2S{S + l))St,Sj)) 



2 , 



+2((>S'^ S^Sl^; SJ)) 2{{Sl{S^ Sj + SjS^); SJ)) 
-B-Gp^'^-^B^GtJ^ + B^Gt;'^ 

-{K,,, + K,,)Glp^ + (2K2,i + AK^,)Gf^'^ - 6K,,,Gf^^^ + AK^.G^ 

uG-/'^ = A7^^^^ + Wj,,{{{S^{6S^S^ -2SiS +l)y,SJ)) 

^ k 

+2((S'j S^ S^; SJ)) — 2{{S^{S- Sj^ + S^ SJ); SJ)) 
+ 2 ~ 

+(7^2,^ + i^4,.)G',7'^ - (27^2,^ + AK^^,)Gp'^"^^ + 67^^,4^,:^.^^^''^ - AK^^.Gp' 
uG^J^'^ = Ajj'^'^ — ^ Jik(^{{{S- Sj + SjSl)Sj; SJ)) — 2{{Sj Sj S^; SJ)) 

k 

~B+Glp^ + 25"G'5+'^ - 47^2,^(^5+'^ - G^/^'^) 
-ir,,4(8G;f'++'^ - 246-1^'++'^ + 326-^++'^ - 16^++'^), 
ujGj^j = + ^ Jjfc(^((S';j (Sf + S^ S^); SJ)) — 2{{S,i S,i SI; SJ 

k 

+B-G;J^^ - 2B-'G;r^^ + AK^^^iG^f'^ - G;r'^^) 



-K,48G-j-^'^''^ - 2AG--^'^''^ + 32G,7"'^ - IQG^f'^), 
uG^'^ = A^J^ + J2jM{{{S^Sr + SrSnSj;Sj)) -3{{Sj:{S^Sj + StSn;Sj)) 



These equations are exact. The important point now is that the anisotropy terms in 
these equations can be simplified by using formulae which reduce products of spin 



operators by one order. Such relations were derived in Ref. [14 



2S-m 



2S-m 
i=0 

AS,m) 



The coefficients 6l ' are tabulated in Ref. |T^ for general spin. For spin 5=1, 



only the coefficients with m = 0,1,2 occur: 6q^''^^ = 62''^^ = 0;6i'''^^ = 1,^0^'^'' 



0,5f'^) = l,5S^'^) = l. 

Application of these relations, effects the reduction of the relevant Green's func- 
tions coming from the anisotropy terms in equations (0): 

















_ ^-{zf,T _ 
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(8) 



The higher Green's functions coming from the anisotropy terms have thus been 
expressed in terms of the lower-order functions already present in the hierarchy; i.e. 
with respect to the anisotropy terms, a closed system of equations of motion results, 
so that no decoupling of these terms is necessary. In other words, the anisotropy is 
treated exactly. For higher spins, S > 1, one can proceed analogously. For this, one 
needs even higher-order Green's functions but again, applying equations (|^ reduces 
the relevant Green's functions by one order, which in turn leads to a closed system 
of equations obviating the decoupling of terms coming from the anisotropics. 

No such procedure is available for the exchange interaction terms, however, so 
that these still have to be decoupled. For spin S' = 1, we use RPA-like approxima- 
tions to effect the decoupling: 

{{Srs',;Sj)) ^ {Sr)Gtf + {S',)Gtf 
{{S^S^Sr,SJ)) ^ {S^)G^/'^ + {S^S^)Gtf. (9) 

Note that we have constructed the decoupling so as not to break correlations having 
equal indices, since the corresponding operators form the algebra characterizing the 
group for a spin S = 1 system. For spin 5 = 1, this decoupling model leads to 8 
diagonal correlations for each layer i: 

These have to be determined by 8 x equations, where A^ is the number of layers. 
We have not attempted to go beyond the RPA- approximation because a previous 
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comparison of Green's function theory with 'exact' quantum Monte Carlo calcula- 
tions for a Heisenberg hamiltonian for a monolayer with S = 1/2 in a magnetic field 
showed RPA to be a remarkably good approximation 

We now apply the above reduction, Eqs. (||), and the decoupling of the exchange 
interaction terms, Eqs. (^, to the monolayer with spin S = 1. Then, after per- 
forming a two-dimensional Fourier transformation, one obtains a set of equations of 
motion, which, in compact matrix notation (dropping the layer index), is as follows: 

(tul -r)G^ = A^, (10) 

where G"^ and are 8-dimensional vectors with components G"'^ and A°'^ where 

a = +,—, z, z+, —z, ++, , zz, and 1 is the unit matrix. The 8x8 non-symmetric 

matrix T is given by 



r = 



/ 


Hi 





-Ht 


K2 














\ 







-Hi 


Hk 





-K2 















-^Hk 


























K2 - ^{6S^S^ -A) 


-{S+S+)Jk 


{{2S- - l)S+).h 







-H' 











{s-s-)Jk 


-K2 + ^(6S^S^ -4) 


-(5-(2S^-l))Jfc 





-H^ 





H+ 


^2H~ 






-{{2S^ -l)S+)Jk 





2{S+S+)Jk 


-H+ 





2H^ 















(5-(25^ -l)>Jfc 


-2(S-S-)Jk 





H~ 





-2H^ 







V 


3(5-(25^-l))Jfc 


-3((2S^ -l)S+>Jfc 





~3H- 


3H+ 







(11) 





/ 



with the abbreviations 

= 5- + (5")J(g-7k), a = +,-,z 
^ B'' + {S'')Jq, a = +,-,z 
Jk = Jl^., (12) 
K2 = K2 + Ki. 

For a square lattice with a lattice constant taken to be unity, 7k = 2 (cos A;^. + cos ky), 
and g = 4, the number of nearest neighbours. For spin 5* = 1 and S = 3/2, 
the Ki term in the Hamiltonian leads only to a renormalization of the second- 
order anisotropy coefficient: K2{S = 1) = + ^4, and K2{S = 3/2) = K2 + 
respectively. Only in the case of higher spins, S" > 2, are there non-trivial corrections 
due to the fourth- order anisotropy coefficient. 

If the theory is formulated only in terms of G^, there is no equation for determin- 
ing the {S^S^) occuring in the F— matrix. It is for this reason that we introduced 

in Eq. (^), for which the F— matrix turns out to be the same, so that, in general, 
one can take a linear combination of G"*" and G~ and their corresponding inhomo- 
geneities: 

G = (l-a)G- + aG+, 
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A = (l-a)A~ + aA+. 



(13) 



Hence, the equations of motion are 

{col - r)G 



(14) 



c 



from which the desired correlations C = (1 — a)C~ + aC~^ can be determined. The 
parameter a is arbitrary (0 < a < 1). The correlation vector for spin 5" = 1 in terms 
of the 8 correlations mentioned above is 

/ {{{l-a)S- + aS+)S+) ^ 

{{{l-a)S- + aS+)S-) 
(((l-a)5- + a5+)5") 
(((1 - a)S- + aS+){2S'S+ - S+)) 
(((1 - a)S- + aS+){2S-S' - S-)) 
{{{l-a)S- + aS+)S+S+) 
{{{l-a)S- + aS+)S-S-) 
{{{l-a)S- + aS+){6S''S' -4)) J 

(1 - a)(2 - {S') - (5^5^)) + a{S+S+) 
(1 - a){S-S-) + a{2 + {S') - {S'S')) 

{1 - a){S-S') + a{{S'S+) - {S+)) 
(1 - a)(2 + {S') - 3{S'S')) - a{S+S+) 
(1 - a){S-S-) + a{{S') + 3{S'S') - 2) 
(1 - a){2{S+) - 2{S'S+)) + a{S+S+S+) 

(1 - a){S-S-S-) + a2{S-S^) 
a){6{S-S') - 4:{S-)) + a{2{S+) - 6{S'S+)) J 

where one can introduce the identity (for spin S = 1): {S~^S~^S~^) = {S~S~S~) = 0. 
The inhomogeneity vectors for spin S = 1 are given by 



(15) 



( ] 




' 2{S') ^ 




^ A+'+ ^ 




( \ 











A-^+ 




-2{S') 


A^- 








A''+ 










6{S'S') - 4 




A'+'+ 




2{S+S+) 


j^-z,- 




-2{S-S-) 




j^-z,+ 




4 - 6{S'S^) 






A{S'S+) - 2{S+) 











'- 









A— '+ 




2{S-) - A{S-S') 






^ 6{S-) - 12{S-S^) ^ 








^ 12{S'S+) - 6{S+) J 



(16) 
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The correlations C are related to the Green's functions via the spectral theorem. In 
order to determine these, we apply the eigenvector method already used in Ref. 
and explained there in detail. This method is quite general and not restricted to 
the 8x8 problem above; it also makes the extension of the theory to multi-layer 
systems tractable. 

The essential steps in deriving the coupled integral equations for determining 
the correlations C are now outlined. One starts by diagonalizing the non- symmetric 
matrix T of equation (|Tip 

LTR = n, (17) 

where R is a matrix whose columns are the right eigenvectors of T and its inverse 
L = R'-*- contains the left eigenvectors as rows, where RL = LR = 1. Multiplying 
Eq. (|T4|) from the left by L and inserting RL = 1 yields 

{uji - n)g = A, (18) 

where we introduce Q = LG and A = LA. Here ^ is a new vector of Green's 
functions with the property that each component Q'^ has but a single pole 

A 



UJ — UJt 



(19) 



This allows the application of the spectral theorem [0 to each component separately, 
with C = LC: 

r = 4^ + T^^- (20) 

V = LD is the correction to the spectral theorem needed in case there are vanishing 
eigenvalues. The corresponding components are obtained from the anticommu- 
tator Green's function Qrj=+i- 



= \ hm 0.^,=+! = \ hm ^^^^ = b^^oAUi , (21) 



i.e. T>'^ is non-zero only for eigenvalues lOr = 0. Denoting these by and the 
corresponding left eigenvectors by L°, one obtains from the Eq. (|21|) 



1^' = l-K=+i = ^L°A,=+i = iL°(A,=_i + 2C) = C°. (22) 

Here, we have exploited the fact that the commutator Green's function is regular at 
the origin (called the regularity condition in [^): 

L°A,=-i=E^°^:;=-i = 0. (23) 

a 
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The desired correlation vector C is now obtained by multiplying the correlation 
vector C, Eq. (]20|), from the left by R: 



C = R £: L A + R"L"C. (24) 

Here, the two terms on the right-hand side belong to the non-zero and zero eigen- 
values of the F— matrix, respectively. R is the matrix whose columns are the right 
eigenvectors of the F-matrix with eigenvalues Ur ^ and L is the corresponding 
matrix whose rows are the left eigenvectors with eigenvalues uJr ^ 0. £^ is a diagonal 
matrix whose elements are the functions — . , . The matrices HP and L° con- 

cxp (pojt) — 1 

sist of the right (column) and left (row) eigenvectors corresponding to eigenvalues 
Ur = 0. This constitutes a system of integral equations which has to be solved 
self-consistently. 

Note that the right-hand side of Eq. (p^) contains a Fourier transformation, 
which can be made manifest by writing out the equations for each component of C 
explicitly: 

-I .jj- ..jY N n n m 

C^ = - dkj dkyY,{T.llR^3^3k5,kLuAi + Y.R%LlC,). (25) 
^ •'^ 1=1 j=ik=i j=i 

Here we have i = 1,..,N correlations Cj corresponding to the N-dimensional F- 
matrix with n non-zero and m zero eigenvalues {n + m = N). The momentum 
integral goes over the first Brillouin zone. For the case of a monolayer with spin 
S = 1, the total number of eigenvalues is = 8, and one can show, by writing down 
the characteristic equation of the F— matrix, that 2 eigenvalues are exactly zero; i.e. 
n = 6, m = 2. 

In general this matrix equation can be ill-defined, for, without loss of gener- 
ality, one can choose the field component to be zero, in which case the cor- 
relations {{S^)"^{S^)"') are the same as {{S~)"'{S^)"^) . This leads to a system of 
overdetermined equations. These equations are solved by means of a singular value 
decomposition [T^, which is now illustrated for spin S = 1. In this case, we have 
{S+) = {S-), {S+S+) = {S-S-), and {S^S+) = {S-S'); i.e. there are only 5 inde- 
pendent variables defining the 8 correlations C. We denote these variables by the 
vector V 

/ is-) \ 

{S-S-) (26) 
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Then, the correlations C can be expressed as 

C = U° + UcV 



(27) 



with 



and 



V 






—a 



2 - 2a 


6a -4 



^ 2 - 2a \ 
2a 


2 - 2a 
-2a 






a- 1 
a 


1-a 
a 






a 

1-a 


—a 
1 — a 








1 



2a -2 

2a 
6 - 12a 



a- 1 
—a 


3a — 3 
3a 






(2^ 



(29) 



Now we write the 8x5 matrix Uc in terms of its singular value decomposition: 



Up = UwV, 



(30) 



where w is a 5 x 5 diagonal matrix whose elements are referred to as the singular 
values. These are in general zero or positive but in our case they are all > for 
< a < 1. U is an orthogonal 8x5 matrix and V is a 5 x 5 orthogonal matrix. 
From Eqs. (|2|) and (|7|) it follows that 



UeV = R ^ L A + R°L°(UeV + u°) - u°. 



(31) 



To get V from this equation, we need only multiply through by u^. = Vw """U, 
which yields the system of coupled integral equations 



V = ur^fRf L A + R°L°(ueV + u°) - u° 



(32) 
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or more explicitly with i = 1 , . . . , 5 




2 5 8 

+ T.RlMJ2i^c),pv, + iu%} - Y.K'Mu%. (33) 

1=1 p=l k=l 

This set of equations is not overdetermined (5 equations for 5 unknowns in the 
present example ) and is solved by the curve-following method described in Appendix 
A. 



3 The multilayer case 



For a ferromagnetic film with layers and spin S = 1 one obtains 8A^ equations of 
motion for the 8A^-dimensional Green's function vector G 



(34) 



where 1 is the 8A^ x 8A^ unit matrix, and the Green's function and inhomogene- 
ity vectors consist of A^ 8-dimensional subvectors which are characterized by layer 
indices i and j 



1 



A'^^ik.u) = {l-a^A^J^ + aA^f. (35) 
The equations of motion are then expressed in terms of these layer vectors, and 8x8 



submatrices Tij of the 8A^ x 8A^ matrix F 



( Tn 

r2i 



ri2 

^22 



'^2N 



\ 



Nl 



N2 



N 



(36) 

In the multilayer case, the F matrix reduces to a band matrix with zeros in the Fij 
sub-matrices, when j > i + 1 and j < i — The diagonal sub-matrices Fa are of size 
8x8 and have the same structure as the matrix which characterizes the monolayer, 
see Eq. (|n|). The matrix elements of Tu contain terms depending on the layer index 
i and additional terms due to the exchange interaction between the atomic layers. 



Tja 

m 



Bt + {St)Uq - 7k) + J^,^+l{S?+l) + J^,-l{St^ 



a 



a 



+, 



(37) 
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The dipole coupling is treated in the mean field limit, which was shown to be a good 
approximation for coupling strengths much weaker than the exchange coupling^]. 
In this case, the dipole terms make additive contributions to the magnetic field 
components i?f , 

N 



N 



where the lattice sums for a two-dimensional square lattice are given by 



(38) 



(39) 



where n = \i — j\. The indices {Im) run over all sites of the square jth layer, 
excluding the terms with + + = 0. For the monolayer (A^ = 1), one has 
i = j, and obtains in particular T° ~ 4.5165. As can be seen from Eqs. (PBD, the 
dipole coupling reduces the effect of the external field component in 2;-direction and 
enhances the effect of the transverse field components B^. 

The 8x8 non-diagonal sub-matrices Tij for j = i ±1 are of the form 










(St) 

















\ 







(Sf> 


-{SD 















































-i(6Sf5f-4> 




((2Sf - 























+ i(6Sf5f -4) 


(-5-(2Sf -1)> 





















-({25f - 1)5+) 
























V 





{S-i2S!-l)) 


-2{S-S-) 



















3(S-(2Sf -1)> 


-3((2Sf - l)S+> 




















/ 



(40) 



We now demonstrate that, if there is an eigenvector L° with eigenvalue zero 
for the sub-matrix Fii, then there is also a left eigenvector of F corresponding to 
eigenvalue zero with the structure 



(0, 



0,L^,0, 



.0) 



where, for spin 5* 



(rO rO rO rO rO rO rO rO 

1-*^?.,+ ) ^i,-^ ^i,z^ ^i,z+J ^i,-zJ ^i, ' ^i,^ 



(41) 



(42) 



Multiplying F from the left by L° results in products of L,^ with sub- matrices Fj^. 
The product with Fa must be zero, since the diagonal blocks of F have the same 
structure as the monolayer matrix, Eq. (|n|). For the off-diagonal blocks, Fjj, the 
product is also zero because of the regularity conditions for layer i, derived from 
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the fact that the commutator Green's functions have to be regular at the origin; see 
Refs. P, § : 

L^ia = 

E Ll Ar = 0, (43) 

a 

E Li Ar = 0. 

a 

Multiplying the non-diagonal matrix (|40|) from the left by the eigenvector (|4^) and 
then applying the regularity conditions Eqs. (0) yields zero. Hence, we have shown 
that there are as many zero eigenvalues of F as there are zero eigenvalues of all of 
the diagaonal blocks Fa. Since each diagonal block has 2 zero eigenvalues (because 
each block has the same structure as the monolayer matrix), there must be at least 
2N zero eigenvalues of the matrix T. 

Therefore, apart from dimension, the equations determining the correlation func- 
tions for the multi-layer system have the same form as for the monolayer case: 

C = R^LA + R°L°C. (44) 

The matrices R and L have to be constructed from the right and left eigenvectors 
corresponding to non-zero eigenvalues as before, whereas the matrices R° and L° 
are constructed from the eigenvectors with eigenvalues zero. 



4 Numerical results 

The results of the numerical calculations presented in this paper are meant to demon- 
strate that our formulation in handling the single-ion anisotropy works in practice. 
To this end we take the magnetic field components and the dipole coupling constant 
to be zero and investigate the magnetization as a function of the anisotropy strength 
and the temperature for a square monolayer with spin S = 1. In this case there is 
only a magnetization (S^) in z-direction. 

In Fig. 1 we show results of mean field (MFT) calculations for (S^) and (S^S^) as 
a function of the temperature for different anisotropics in the range of < < 300 
obtained in two ways. The first is an exact diagonalization of the mean field Hamil- 
tonian, which is possible because of its one-body nature. If our Green's function 
theory (GET) for the anisotropy term is exact, calculations with the Green's function 
program in the mean field limit (no momentum dependence on the lattice: 7k = 
of Eq. (p!2D) should give identical results. This is indeed the case; both results are 
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S=1 monolayer: mean field theory 




temperature 

Figure 1: Results of mean field calculations using either the mean field limit (7k = 0) 
of the Green's function program or an exact diagonalization of the corresponding mean 
field Hamiltonian. Both results are identical. For a monolayer with S = 1, the magneti- 
zation component (S^) and the correlation (S^S^) in MFT are shown as functions of the 
temperature for anisotropy coefficients in the range < K2 < 300; the exchange coupling 
strength is J = 100. 



MFT: Anderson Callen decoupling 
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Figure 2: The figure displays results of the MFT limit of a GFT with the Anderson- 
Callen decoupling, demonstrating the difference from the exact results of Fig. 1. The 
magnetization (S^) and correlation {S^S^) arc shown only up to K2 = 100, where already 
the differences arc large; results for K2 > 100 lie outside the temperature scale of the 
figure. Note that the values for (S^S^) = 2/3 contrast with the exact results. 
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indistinguishable in Fig. 1. The precise agreement of these very different methods 
of calculations provides a check on the numerical procedures. 

Fig. 2 presents the results of a MFT calculation using the Anderson-Callen de- 
coupling of the single-ion anisotropy terms. The shortcoming of this decoupling 
is seen by comparing with the exact results of Fig. 1. One observes that, up to 
K2 — 10, the approximate calculation overshoots the exact one only slightly, but 
with increasing K2 the disagreement becomes worse and worse. The results for 
K2 > 100 lie outside the frame of the figure. 

In the MFT results of Figs. 1 and 2 the well-known shortcoming of MFT is 
evident, the violation of the Mermin- Wagner theorem: there is a finite Curie tem- 
perature for vanishing anisotropy: TqI^I^{K2 = 0) = |J = 133.33 for an exchange 
coupling strength of J = 100. For arbitrarily large values of the anisotropy, the Curie 
temperature in MFT is obtained analytically: T^J^^{K2 00) = S'^qJ / {S{S+1)) = 
200 for S = 1, J = 100 and g = 4 (q is the coordination number of a square lattice). 
This limit is almost reached numerically for K2 = 300 as can be seen in Fig. 1. 

Our Creen's function theory with the RPA-like treatment of the exchange terms 
fulfills the Mermin- Wagner theorem: Tcurie ~^ for K2 0. 



GFT: S=1 monolayer 
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Figure 3: Results of the Green's function theory (with RPA-like decouplings of the 
exchange terms) with the same input as in Fig.l are shown for {S^) and (S^S^) as functions 
of the temperature for various anisotropy coefficients K2. Note the significant differences 
from the mean field results of Fig.l. 

Comparison of the MFT results of Fig. 1 with the GFT results in Fig. 3 reveals 
major differences between MFT and GFT with respect to the temperature depen- 
dence of (S^) for different anisotropics K2, particularly in the low temperature region 
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and for small anisotropies. For large anisotropies it can be shown analytically that 
the full Green's function theory approaches the MFT limit, Tcurie = 200, when the 
anisotropy becomes arbitrarily large (see Appendix B). This is physically reasonable 
because, in the large anisotropy limit, GFT approaches the Ising limit, and, for the 
Ising model, a RPA treatment is identical with the mean field approach. The results 
of the exact treatment of the single-ion anisotropy term shown in Fig. 3 represent a 
significant improvement over the decoupling of this term proposed by Anderson and 
Callen |T^ and the different decoupling of Lines|T^, both of which yield a diverging 
Curie temperature Tcurie ~^ oo for K2 00. (See also Appendix B of Ref. in 
this connection.) 



Curie temperatures: MFT and GFT 
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Figure 4: Comparison of the Curie temperatures calculated with the present exact 
treatment of the anisotropy terms, the Anderson-Callen decoupling|^] and MFT. The first 
two approaches fulfill the Mermin- Wagner theorem: Tc for K2 0, whereas the 
MFT result does not. For large anisotropies {K2 ^00), the new model approaches slowly 
the mean field results, as it should do, whereas the Anderson-Callen decoupling procedure 
leads to a diverging Tq 



To show the difference between the new model and the Anderson-Callen decou- 
pling more clearly, we compare in Fig. 4 the Curie temperatures obtained from 
MFT, the new Green's function theory, and the Green's function theory with the 
Anderson-Callen decoupling of Refs. For small anisotropies, there is only 

a slight difference between the two GFT results which, in contrast to MFT, obey 
the Mermin- Wagner theorem. However, on increasing the anisotropies, the GFT 
results deviate from one another significantly: for K2 —>■ 00, the Anderson-Callen 
result diverges, whereas the exact treatment approaches the MFT limit, albeit very 
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slowly. 

The difference between the exact treatment of the anisotropy terms and the 
approximate Anderson-Callen decoupling is further demonstrated in Fig. 5, where 
the magnetizations (S^) as a function of the temperature for different values of K2 
are compared. We see that, for small anisotropics, there is rather good agreement, 
which, however, worsens as K2 increases. Another difference concerns the second 
moments (S^S^), which, in the case of the Anderson-Callen decoupling, approach 
the value {S^S^)iT Tcwie) = 2/3 (see Ref. 0), whereas in the exact treatment, 
the values of {S^S^){T Tcurie) are larger than 2/3. This is as it should be, 
because, as shown in Appendix B, {S^S^) — > 1 for K2 — > 00. 

K2: exact (open circles); decoupling (lines) 
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Figure 5: Comparison of GFT calculations for (S^) and {S^S^) as functions of the 
temperature using the exact treatment of the anisotropy (open circles) and the Anderson- 
Callen decoupling used in in Refs. II, I (small dots). 



5 Conclusions 

We have presented a formal theory for the magnetization of thin ferromagnetic films 
on the basis of many-body Green's function theory within a Heisenberg model with 
anisotropics. The essential improvement over our previous work^, ^ is the exact 
treatment of the single-ion anisotropy brought about by the introduction of higher- 
order Green's functions. Previously, the anisotropy term was treated by approximate 
decoupling procedures only at the level of the lowest-order Green's functions. The 
exchange interaction terms are decoupled using an RPA-like approach. We did not 
try to go beyond RPA since our comparison with 'exact' quantum Monte Carlo 
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results has shown this to be a very good approximation ||l5 



Numerical calculations of the magnetization as a function of the temperature 
for various anisotropics K2 (no external field, no dipole-coupling) demonstrate the 
superiority of the new spin wave model over MFT. In particular, there is no viola- 
tion of the Mermin- Wagner theorem. The Anderson-Callen decoupling used in our 
previous work gives results close to those of the new model when the anisotropy K2 
is small but, as the anisotropy increases, the difference between the two approaches 
becomes larger: the new model approaches the MFT limit as it should do, whereas 
the Curie temperature from the Anderson-Callen decoupling diverges. 

Our new formulation should allow a future investigation of the reorientation 
problem when switching on the magnetic field and/or the dipole coupling. The 
treatment of multi-layer systems and spin S > 1 should be possible. 

We are indebted to A. Ecker and P.J. Jensen for discussions. 



Appendix A: The curve-following procedure 

Consider a set of n coupled equations characterised by m parameters {Pi]i = 
1,2 ... , m} and n variables {Vi] i = 1,2, . . . , n}: 

S,{P[m];V[n]) = 0,i = l,...,n/. (45) 

In our case, the parameters are the temperature, the magnetic field components, 
the dipole coupling strengths, the anisotropy strengths, etc; the variables are the 
spin-correlations. The coupled equations are obtained by defining the Si from the n 
self-consistency equations for the correlations vector v (Eq. (p^)): 

S = v-u-i(R^LA + R°L°(ucV-u°) -u°). 

For fixed parameters P, we look for solutions Si = at localised points, V[n], 
in the n-dimensional space. If now one of the parameters Pk is considered to be 
an additional variable Vo (in this paper, the temperature is taken as the vari- 
able), then the solutions to the coupled equations define curves in the {n + 1)- 
dimensional space V[n + 1]. From here on, we denote the points in this space by 

{Vi] i = 0,1,2, , n}. The curve-following method is a procedure for generating 

these solution-curves point by point from a few closely-spaced points already on a 
curve; i.e. the method generates a new solution-point from the approximate direc- 
tion of the curve in the vicinity of a new approximate point. This is done by an 
iterative procedure described below. If no points on the curve are known, then an 
approximate solution point and an approximate direction must be estimated before 
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applying the iterative procedure to obtain the first point on the curve. A second 
point can then be obtained in the same fashion. If at least two solution-points are 
available, then the new approximate point can be extrapolated from them and the 
approximate direction can be taken as the tangent to the curve at the last point. 

The iterative procedure for finding a better point, V, from an approximate 
point, V°, is now described. One searches for the isolated solution-point in the 
n-dimensional subspace perpendicular to the approximate direction, which we char- 
acterise by a unit vector, u. The functions Si are expanded up to first order in the 
corrections about the approximate point, V°: 

s,{y)^s,{y°) + Y.^ Av„ (46) 

where AVj — Vj—V°. At the solution, the Si are all zero, whereas at the approximate 
point V° the functions have non-zero values, S°; hence, one must solve for the 
corrections AVj for which the left-hand side in the above equation is zero: 

Ay, = -5°;{i = l,2,...,n}. (47) 

These n equations are supplemented by the constraint requiring the correction to 
be perpendicular to the unit direction vector: 

n 

J2u]AVj^0. (48) 

i=o 

This improvement algorithm in the subspace is repeated until each of the Si is 
sufficiently small. In practice we required that J2i {^iY < e, where we took e = 
10~^^. If there is no convergence, the extrapolation step-size used to obtain the 
original V° is halved, a new extrapolated point obtained, and the improvement 
algorithm repeated. 

The curve-following method is quite general and can be applied to any coupled 
equations characterised by differentiable functions. By utilizing the information 
about the solution at neighbouring points, the method is able to find new solutions 
very efficiently, routinely converging after a few iterations once two starting points 
have been found. 

Appendix B: Curie temperature for K2 ^ oo 

We show analytically that the Curie temperature of the Green's function theory 
with the exact treatment of the anisotropy for a square- lattice monolayer with S=l 
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approaches the mean field value when the anisotropy coefficient goes to infinity, 
whereas the Anderson-Callen decoupling leads to a divergence in this limit. 

For the case of a single magnetic direction, the 8x8 problem of Eq. (|TUp reduces 
to a 2 X 2 problem for the Green's functions Gj^j'^ = {{S^-jSj')) and G^^^'^ = (((2S'f — 
l)Sj'; Sj')) . For this special case, it is possible to derive analytical expressions for 
the correlations {S^) and {S^S^): 
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{s-s^) = 2 - {s') - {s's') = ^ r dK r dky{—^ 

TT'^ Jo Jo ^UJ^ — 

-(2(5^)(^- - {S^)Jo) + M6{S^S^) - 4))^^;^] } (49) 

{S-{2S^ - 1)S+) = (S^) - 1(6(5^^^) -4) = ^ FdK rdky{^^ 

Z 7r^ JO Jo ^ uj^ — u 

mS'S^) - 4)(u;+ - {S^Jo) + ^2(2(5^)))-^^^ 
-mS'S^) - A){uj- - {S^)Jo) + K2{2{S^)))^^^j^] } (50) 



with 



= ^(^^)(^o - Jk) ± JkI - 1{Q{S^S^) - 4)K,A - \{S^)'. (51) 



At the Curie temperature, {S^) — >• 0, so that the equation for a;^ becomes 



uj^({S^) = 0) = ±a;° = ±^jK,{K2 - \{Q{S^S^)) - 4)Jk), (52) 

Equation (|49|) then reduces to 

2 - (S'S^) = {6{S'S') - 4)^ r dK r dky^ coth(/5cu72). (53) 

vr^ JO JO Zuj^ 



For large K2, cu^ = K2^1 - (6(5^5^) - 4) Jk/(2ir2) ^ Passing to the limit 
K2 ^ 00, one obtains from Eq. (|^) at Tcurie 

{S'S') = 1. (54) 

Now, expanding Eq. (^OD around {S^) = 0, and comparing the coefficients of {S^) 
of the resulting equation, one has at TQu^te 

^ = V-lo ^^-io aJ^ )coth(^) 

+ (^ - Jo)/3(6(g-g-) - 4) ^^^^r }• (55) 
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Noting that dk^ Jq dkyj]^ = and that cj" ~ K2 for large K2, one obtains from 
Eq. (H) 

1 - coth(^) = -JomS'S^) - 4))p^^-^. (56) 



Again, goint to the hmit K2 — > 00 and using Eq. (|54D , one obtains for the Curie 
temperature 

Tcurie = Jol{Q{S'S') - 4) = Jo = 4J (57) 

This is just the MET result! This is physically reasonable because a large anisotropy 
approaches the Ising limit, and the RPA for the Ising model is identical to its mean 
field treatment. 

This is in contrast to the result of the decoupling procedure. In Appendix B of 
Ref. we have shown that the Anderson-Callen decoupling of the anisotropy term 
leads for a square monolayer to a Curie temperature 

87rJ/3 



- Curie 



ln(l + 3TtU/K2 
which diverges for K2 —>■ 00! 



(5^ 



22 



References 

[1] R.Sellmann, H. Fritzsche, H. Maletta, V. Leiner, and R. Siebrecht, Phys. Rev. 
B 64 (2001) 054418 

[2] S. Putter, H.F. Ding, Y.T. Millev, H.P. Oepen, and J. Kirschner, Phys. Rev. B 
64 (2001) 092409 

[3] A. Moschel, K.D. Usadel, Phys. Rev. B 49 (1994) 12868; P.J. Jensen, K.H. 
Bennemann, in 'Magnetism and Electronic Correlations in Local-Moment Sys- 
tems: Rare-Earth Elements and Compounds', ed. M. Donath, P.A. Dowben, 
and W. Nolting, (World Scientific, 1998), p. 113-141 

[4] P.J. Jensen, K.H. Bennemann, Sohd State Comm. 100 (1996) 585, ibid. 105 

(1998) 577; A. Hucht, K.D. Usadel, Phys. Rev. B 55 (1997) 12309 

[5] P. Probrich, P.J. Jensen, P.J. Kuntz, Eur. Phys. J. B 13 (2000) 477 

[6] P. Probrich, P.J. Jensen, P.J. Kuntz, A. Ecker, Eur. Phys. J. B 18 (2000) 579 

[7] Wenh Guo, L.P. Shi, D.L. Lin, Phys. Rev. B 62 (2000) 14259 

[8] R. P. Erickson, D. L. Mills, Phys. Rev. B 44 (1991) 11825 

[9] C. Timm, P.J. Jensen, Phys. Rev. B 62 (2000) 5634 

[10] A.B. Maclsaac, K. De'BeU, J.P. Whitehead, Phys. Rev. Lett 80 (1998) 616 

[11] T. Herrmann, M. Potthoff, and W. Nolting, Phys. Rev. B 58 (1998) 831 

[12] J.F. Devlin, Phys. Rev. B 4 (1971) 136 

[13] F.B. Anderson, H.B. Callen, Phys. Rev. 136 (1964) A1068 

[14] P.J. Jensen, F. Aguilera-Granja, Phys. Lett. A 269 (2000) 158 

[15] A.Ecker, P. Probrich, P.J. Jensen, P.J. Kuntz, J. Phys.: Condens. Matter 11 

(1999) 1557 

[16] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical 
Recipes, Cambridge University Press, 1989 

[17] W. Gasser, E. Heiner, K. Elk, in 'Greensche Funktionen in der Festkorper- und 
Vielteilchenphysik', Wiley- VHC, Berlin, 2001, Chapter 3.3 

[18] M.E. Lines, Phys. Rev. 156 (1967) 534 



23 



